05: Diversity metrics, microbiome
Eva Kim
12/11/2024
Descriptive microbiome stats and cleaner code for dissertation / publication.
1 Load libraries and directories
library(magrittr)
library(phyloseq)
library(fs)
library(tidyverse)
library(vegan)
library(rstatix)
library(ggpubr)
library(microViz)
library(pals)
library(RColorBrewer)
library(colorBlindness)
library(microbiome)
library(ggrepel)
library(DT)
library(plotly)
library(ape)
library(picante)
library(phytools)
library(DESeq2)
library(parallel)
library(doBy)
library(ragg)
library(cowplot)
library(ggtext)rm(list=ls()) # clear environment
git.dir = file.path("/hpc/group/dmc/Eva/GLP_KO_Microbiome_final") # CHANGE ME
ps.rds <- file.path(git.dir, "ps.GLP1.rds")
ps <- read_rds(ps.rds)
ps## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3304 taxa and 234 samples ]
## sample_data() Sample Data: [ 234 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq() DNAStringSet: [ 3304 reference sequences ]fig.dir = file.path(git.dir, "figures")
map.file = file.path(git.dir, "Metadata_GLP1.csv")
meta.df = read_csv(map.file, show_col_types = FALSE)Note that this phyloseq object is already filtered such that only samples with reads > 1000 are available. See here:
sample_sums(ps) %>% as.data.frame() %>%
dplyr::rename(Reads = ".") %>%
arrange(Reads) The names are hard to interpret because they are the SRA accession numbers. Let’s change the names in the phyloseq object.
phyloseq::sample_names(ps) <- phyloseq::sample_data(ps)$Label1.1 Change metadata to make variables factors
meta.df$Genotype <- factor(meta.df$Genotype, levels = c("WT", "KO"))
meta.df$Sex <- factor(meta.df$Sex, levels = c("Female", "Male"))
meta.df$Diet <- factor(meta.df$Diet, levels = c("Normal_Chow", "High_Fat_Diet"))
meta.df$Intranasal_Treatment <- factor(meta.df$Intranasal_Treatment, levels = c("PBS", "HDM"))
meta.df$Surgery <- factor(meta.df$Surgery, levels = c("None", "Sham", "VSG"))
meta.df$Group <- factor(meta.df$Group, levels = c("NA", "Control", "1", "2"))meta.df %>% dplyr::filter(Label %in% rownames(phyloseq::sample_data(ps))) %>%
mutate(Label_graph = Label) %>%
column_to_rownames("Label") -> meta.ps
ps## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3304 taxa and 234 samples ]
## sample_data() Sample Data: [ 234 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq() DNAStringSet: [ 3304 reference sequences ]sample_data(ps) <- meta.ps
ps## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3304 taxa and 234 samples ]
## sample_data() Sample Data: [ 234 samples by 40 sample variables ]
## tax_table() Taxonomy Table: [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq() DNAStringSet: [ 3304 reference sequences ]Genotype, Sex, Diet,
Intranasal_Treatment, and Surgery are
appropriately factors now.
1.2 Subset phyloseq objects
# Only include samples
subset_samples(ps, Type == "True Sample") -> ps_sam
# remove samples with potential labeling issues
ps_sam %>% subset_samples(row.names(ps_sam@sam_data) != "294" & row.names(ps_sam@sam_data) != "30") -> ps_sam
# Subset by sample type
ps_sam %>% subset_samples(Sample_type == "Lung Tissue") -> ps.LT.asv
ps_sam %>% subset_samples(Sample_type == "Fecal") -> ps.F.asv
ps_sam## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3304 taxa and 221 samples ]
## sample_data() Sample Data: [ 221 samples by 40 sample variables ]
## tax_table() Taxonomy Table: [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq() DNAStringSet: [ 3304 reference sequences ]ps.LT.asv## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3304 taxa and 44 samples ]
## sample_data() Sample Data: [ 44 samples by 40 sample variables ]
## tax_table() Taxonomy Table: [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq() DNAStringSet: [ 3304 reference sequences ]ps.F.asv## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3304 taxa and 177 samples ]
## sample_data() Sample Data: [ 177 samples by 40 sample variables ]
## tax_table() Taxonomy Table: [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq() DNAStringSet: [ 3304 reference sequences ]# Relative abundance
ps_sam.ts <- transform_sample_counts(ps_sam, function(x) x/sum(x))
ps.LT.asv.ts <- transform_sample_counts(ps.LT.asv, function(x) x/sum(x))
ps.F.asv.ts <- transform_sample_counts(ps.F.asv, function(x) x/sum(x))1.3 Clean-up Metadata
Let’s subset the metadata to those that are in the
ps_sam object.
meta.df %>%
filter(meta.df$Label %in% rownames(ps_sam@sam_data)) -> meta.sam
meta.sam %>% filter(is.na(Genotype) == TRUE & is.na(Intranasal_Treatment) == TRUE & is.na(Surgery) == TRUE)2 Rarefaction curves
rarecurve(as.data.frame(otu_table(ps_sam)), step = 100, cex = 0.5, tidy = TRUE) -> rare_ps #put in ggplot2 friendly dataframe## Warning in rarecurve(as.data.frame(otu_table(ps_sam)), step = 100, cex = 0.5, :
## most observed count data have counts 1, but smallest count is 2meta.sam %>% dplyr::select(Label, Sample_type, Mouse, Week, DNA_ng_ul_raw) %>%
dplyr::rename(Site = Label) %>%
right_join(rare_ps, by = "Site") -> rare_ps
ggplot(data=rare_ps, aes(x=Sample, y=Species, group=Site)) +
geom_line() + theme_bw() + labs(x = "Number of Reads", y= "Species Richness") +
facet_wrap(~Sample_type, scales = "free") -> p
p
ggplotly(p, tooltip = "Site")3 Relative abundance plots
ps_sam %>% tax_glom("Phylum") %>%
transform_sample_counts(function(x) x/sum(x)) -> ps.phylum.ts
taxa_names(ps.phylum.ts) <- tax_table(ps.phylum.ts)[, 'Phylum']
ps.phylum.ts## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 20 taxa and 221 samples ]
## sample_data() Sample Data: [ 221 samples by 40 sample variables ]
## tax_table() Taxonomy Table: [ 20 taxa by 6 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 20 tips and 19 internal nodes ]
## refseq() DNAStringSet: [ 20 reference sequences ]plot_bar(ps.phylum.ts,x = "Label_graph", fill="Phylum") +
geom_bar(aes(color=Phylum, fill=Phylum), stat="identity", position="stack") +
theme_bw() +
labs(y="Relative Abundance", x = "Sample") +
coord_cartesian(ylim = c(0,1), expand=0) +
scale_color_manual(values=SteppedSequential5Steps) +
scale_fill_manual(values=SteppedSequential5Steps) +
theme(axis.text.x=element_blank()) +
facet_wrap(~Sample_type, scales = "free") -> p
p
ps_melt(ps.phylum.ts) -> phylum.ts.melt
phylum.ts.melt %>%
filter(Sample_type == "Lung Tissue") %>%
doBy::summary_by(Abundance ~ OTU, FUN = median) %>%
arrange(desc(Abundance.median))phylum.ts.melt %>%
filter(Sample_type == "Fecal") %>%
doBy::summary_by(Abundance ~ OTU, FUN = median) %>%
arrange(desc(Abundance.median))ps.phylum.ts %>%
subset_samples(Sample_type == "Fecal") %>%
ps_arrange(desc(Firmicutes), desc(Bacteroidota), .target = "otu_table") -> ps.fecal.phylum.ts
sample_data(ps.fecal.phylum.ts)$Sample_Order <- c(1:nrow(sample_data(ps.fecal.phylum.ts)))
ps.fecal.phylum.ts## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 20 taxa and 177 samples ]
## sample_data() Sample Data: [ 177 samples by 41 sample variables ]
## tax_table() Taxonomy Table: [ 20 taxa by 6 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 20 tips and 19 internal nodes ]
## refseq() DNAStringSet: [ 20 reference sequences ]ps.phylum.ts %>%
subset_samples(Sample_type == "Lung Tissue") %>%
ps_arrange(desc(Proteobacteria), desc(Actinobacteriota), .target = "otu_table") -> ps.lt.phylum.ts
sample_data(ps.lt.phylum.ts)$Sample_Order <- c(1:nrow(sample_data(ps.lt.phylum.ts)))
ps.lt.phylum.ts## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 20 taxa and 44 samples ]
## sample_data() Sample Data: [ 44 samples by 41 sample variables ]
## tax_table() Taxonomy Table: [ 20 taxa by 6 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 20 tips and 19 internal nodes ]
## refseq() DNAStringSet: [ 20 reference sequences ]ps_melt(ps.fecal.phylum.ts) %>%
ggplot(aes(fill=Phylum, y=Abundance, x=Sample_Order)) +
geom_bar(position="fill", stat="identity", width = 1) +
theme_bw() +
labs(y="Relative Abundance", x = "Sample") +
coord_cartesian(ylim = c(0,1), expand=0) +
scale_color_manual(values=SteppedSequential5Steps) +
scale_fill_manual(values=SteppedSequential5Steps) +
facet_wrap(~Sample_type) +
theme(axis.text.x=element_blank()) -> p.fecal
p.fecal
ps_melt(ps.lt.phylum.ts) %>%
ggplot(aes(fill=Phylum, y=Abundance, x=Sample_Order)) +
geom_bar(position="fill", stat="identity", width = 1) +
theme_bw() +
labs(y="Relative Abundance", x = "Sample") +
coord_cartesian(ylim = c(0,1), expand=0) +
scale_color_manual(values=SteppedSequential5Steps) +
scale_fill_manual(values=SteppedSequential5Steps) +
facet_wrap(~Sample_type) +
theme(axis.text.x=element_blank()) -> p.lt
p.lt
ragg::agg_jpeg(file.path(fig.dir, "fig_RelAb_fecal.jpeg"), width = 10, height = 6, units = "in", res = 600)
p.fecal
dev.off()## png
## 2ragg::agg_jpeg(file.path(fig.dir, "fig_RelAb_lt.jpeg"), width = 5, height = 6, units = "in", res = 600)
p.lt
dev.off()## png
## 24 Alpha Diversity: LT Samples
4.1 Calculate
ps.LT.asv %>% estimate_richness(measures = c("Observed", "Shannon")) -> LT_richness
LT_richness$Faith <- pd(otu_table(ps.LT.asv), phy_tree(ps.LT.asv), include.root = FALSE)$PD
LT_richness %>% rownames_to_column(var = "Label") -> LT_richness
gsub("\\.", "-", LT_richness$Label) -> LT_richness$Label
meta.df %>%
filter(meta.df$Label %in% LT_richness$Label) %>% # merge with metadata
select(c(Label:DNA_ng_ul_raw, "Group")) %>%
right_join(LT_richness, by = "Label") -> LT_richness
LT_richness %>% pivot_longer(
cols = contains(c("Observed", "Shannon", "Faith")),
names_to = c("Measurement")
) -> LT_richness_stat # dataframe that works well for ggplot2() and anova_test()
LT_richness_stat$Mouse <- factor(LT_richness_stat$Mouse)
LT_richness_stat$Measurement <- factor(LT_richness_stat$Measurement)
LT_richness_stat$Cohort <- factor(LT_richness_stat$Cohort)LT_richness %>% nrow()## [1] 44LT_richness_stat %>% nrow()## [1] 13244*3## [1] 1324.2 No surgery mice
We are most interested in intranasal treatment, diet, and genotype if possible. We also note that there were some mice whose lung tissues were collected on week 10 vs. week 13.
LT_richness %>% filter(Surgery == "None") -> lt.richness.nosurgery
lt.richness.nosurgery %>%
dplyr::count(Week, Genotype, Intranasal_Treatment) %>% arrange(n)lt.richness.nosurgery %>%
dplyr::count(Week, Intranasal_Treatment, Diet) %>% arrange(n)Let’s graph:
LT_richness_stat %>% filter(Surgery == "None") -> lt.richness.nosurgery.long
lt.richness.nosurgery.long %>%
ggplot(aes(x = Intranasal_Treatment, y = value)) +
geom_boxplot() + geom_point() +
facet_wrap(vars(Measurement), scales = "free") + theme_bw()
lt.richness.nosurgery.long %>%
ggplot(aes(x = Genotype, y = value)) +
geom_boxplot() + geom_point() +
facet_wrap(vars(Measurement), scales = "free") + theme_bw()
lt.richness.nosurgery.long %>%
ggplot(aes(x = Diet, y = value)) +
geom_boxplot() + geom_point() +
facet_wrap(vars(Measurement), scales = "free") + theme_bw()
lt.richness.nosurgery.long %>%
ggplot(aes(x = as.character(Week), y = value)) +
geom_boxplot() + geom_point() +
facet_wrap(vars(Measurement), scales = "free") + theme_bw()
lt.richness.nosurgery.long %>%
ggplot(aes(x = Genotype, y = value)) +
geom_boxplot() + geom_point() +
facet_grid(vars(Measurement), vars(Intranasal_Treatment),
scales = "free") + theme_bw()
lt.richness.nosurgery.long %>%
ggplot(aes(x = Diet, y = value)) +
geom_boxplot() + geom_point() +
facet_grid(vars(Measurement), vars(Intranasal_Treatment),
scales = "free") + theme_bw()
lt.richness.nosurgery.long %>%
ggplot(aes(x = as.character(Week), y = value)) +
geom_boxplot() + geom_point() +
facet_grid(vars(Measurement), vars(Intranasal_Treatment),
scales = "free") + theme_bw()
Graphically, it looks like intranasal treatment, week, and genotype may be of interest.
lt.richness.nosurgery %>% nrow()## [1] 26lt.richness.nosurgery %>%
dplyr::count(Week, Intranasal_Treatment) %>% arrange(n)lt.richness.nosurgery.long %>%
group_by(Week, Intranasal_Treatment, Measurement) %>%
shapiro_test(value) %>%
add_significance()Mainly indistinguishable from normal distribution.
lt.richness.nosurgery.long %>%
group_by(Week, Intranasal_Treatment, Measurement) %>%
identify_outliers(value) -> lt.nosurgery.outliers
lt.nosurgery.outliersOutliers only in week 10 groups.
meta.sam %>% filter(Sample_type == "Lung Tissue") %>%
filter(Label %in% lt.nosurgery.outliers$Label) %>%
dplyr::select(Label, Sample_type:DNA_ng_ul_raw, Notes)Interesting, maybe genotype or sex is playing a role here?
lt.richness.nosurgery.long %>%
ggqqplot(x = "value") + facet_grid(vars(Measurement), vars(Intranasal_Treatment), scales = "free")
lt.richness.nosurgery.long %>%
group_by(Measurement) %>%
anova_test(
dv = value,
between = c(Week, Intranasal_Treatment, Genotype),
) %>% get_anova_table() %>%
adjust_pvalue(method = "holm") %>%
add_significance()Nothing is significant, sadly. Let’s graph with intranasal treatment, as that had the greatest effect on host respiratory health metrics.
lt.richness.nosurgery.long %>%
group_by(Measurement) %>%
t_test(value ~ Intranasal_Treatment) %>%
adjust_pvalue(method = "holm") %>%
add_significance() -> lt.richness.nosurgery.res
lt.richness.nosurgery.res# Change facet label names
alpha.labs <- c("Faith's PD", "Observed Richness", "Shannon")
names(alpha.labs) <- c("Faith", "Observed", "Shannon")
#lt.richness.nosurgery.res$y.position <- c(40, 225, 4.7)
lt.richness.nosurgery.long %>%
ggplot(aes(x = Intranasal_Treatment, y = value, color = Intranasal_Treatment))+
geom_boxplot() +
geom_point() +
facet_wrap(~Measurement, scales = "free_y",
labeller = labeller(Measurement = alpha.labs)) + expand_limits(y = 0) +
theme_bw(base_size = 14) + scale_fill_brewer(palette="RdBu") +
labs(y="Alpha Diversity Value", x="Intranasal Treatment", color = "Intranasal Treatment") +
theme(legend.position = "top") -> lt.alpha.group1.week
# stat_pvalue_manual(lt.richness.nosurgery.res, label = "p.adj.signif") -> lt.alpha.group1.week
lt.alpha.group1.week
4.3 HFD: effect of surgery?
For the after surgery timepoint for lung tissue, look at week 13 only and HFD.
LT_richness %>% filter(Diet == "High_Fat_Diet" & Week == 13) %>%
dplyr::count(Surgery, Intranasal_Treatment)LT_richness %>% filter(Diet == "High_Fat_Diet" & Week == 13) %>%
dplyr::count(Surgery)LT_richness %>% filter(Diet == "High_Fat_Diet" & Week == 13) %>%
dplyr::count(Intranasal_Treatment)LT_richness %>% filter(Diet == "High_Fat_Diet" & Week == 13)%>% nrow()## [1] 23LT_richness_stat %>% filter(Diet == "High_Fat_Diet" & Week == 13) -> lt.richness.hfd
lt.richness.hfd %>%
group_by(Measurement, Surgery) %>%
shapiro_test(value) %>%
add_significance()lt.richness.hfd %>% ggqqplot(x = "value") + facet_grid(vars(Measurement), vars(Surgery), scales = "free")
It looks like this is due to an outlier, but rarefaction curves looked like they plateaued. Proceed with caution.
lt.richness.hfd %>%
group_by(Measurement, Surgery) %>%
shapiro_test(value) %>%
add_significance()lt.richness.hfd %>%
group_by(Measurement, Surgery) %>%
identify_outliers(value) lt.richness.hfd %>%
ggplot(aes(x = Surgery, y = value)) +
geom_boxplot() + geom_point() +
facet_wrap(vars(Measurement), scales = "free") + theme_bw()
lt.richness.hfd %>%
group_by(Measurement) %>%
anova_test(
dv = value,
between = c(Surgery, Intranasal_Treatment),
) %>% get_anova_table() %>%
adjust_pvalue(method = "holm") %>%
add_significance()Nothing is significant. Let’s make a graph showing surgery:
lt.richness.hfd %>%
ggplot(aes(x = Surgery, y = value, color = Surgery)) +
geom_boxplot() + geom_point() +
facet_wrap(~Measurement, scales = "free_y",
labeller = labeller(Measurement = alpha.labs)) + expand_limits(y = 0) +
theme_bw(base_size = 14) + scale_color_brewer(palette="Dark2") +
labs(y="Alpha Diversity Value", x="Surgery", color = "Surgery") +
theme(legend.position = "top") -> fig.lt.hfd
fig.lt.hfd
4.4 Summary figure
plot_grid(
lt.alpha.group1.week, fig.lt.hfd,
labels = "AUTO", ncol = 2, label_size = 20, scale = 0.95,
rel_widths = c(1, 1.3)
) -> fig.lt.alpha
fig.lt.alpha
ragg::agg_jpeg(file.path(fig.dir, "fig_lt_alpha.jpeg"), width = 14, height = 5, units = "in", res = 600, scaling = 1)
fig.lt.alpha
dev.off()## png
## 25 Lung tissue beta diversity
5.1 No surgery mice
For PERMANOVA, let’s have n>=5 for each group.
ps.LT.asv.ts %>% subset_samples(Surgery == "None") -> ps.lt.ts.nosurgery
# Calculate ditances
ps.lt.ts.nosurgery %>% phyloseq::distance(method = "unifrac") -> LT_nosurgery_uni
ps.lt.ts.nosurgery %>% phyloseq::distance(method = "wunifrac") -> LT_nosurgery_weighted_uni
# Ordinate distances
LT_nosurgery.uni_ord <- ordinate(ps.lt.ts.nosurgery, method = "PCoA", distance = "unifrac")
LT_nosurgery.weighted_uni_ord <- ordinate(ps.lt.ts.nosurgery, method = "PCoA", distance = "wunifrac")
# Extract dataframe
ps.lt.ts.nosurgery@sam_data %>% data.frame() -> LT_nosurgery.df
LT_nosurgery.df %>% dplyr::count(Genotype, Intranasal_Treatment) %>% arrange(n)# Set seed for reproducibility
# UniFrac
set.seed(123)
adonis2(formula = LT_nosurgery_uni ~ Intranasal_Treatment * Genotype, data = LT_nosurgery.df, permutations = 9999, parallel = getOption("mc.cores"))# Weighted UniFrac
set.seed(123)
adonis2(formula = LT_nosurgery_weighted_uni ~ Intranasal_Treatment * Genotype, data = LT_nosurgery.df, permutations = 9999, parallel = getOption("mc.cores")) Not significant.
5.1.1 Graphs
plot_ordination(ps.lt.ts.nosurgery, LT_nosurgery.uni_ord, type = "samples", color = "Intranasal_Treatment") + theme_bw() +
labs(title = "UniFrac") +
scale_color_brewer(palette="Dark2")
plot_ordination(ps.lt.ts.nosurgery, LT_nosurgery.weighted_uni_ord, type = "samples", color = "Intranasal_Treatment") + theme_bw() +
labs(title = "Weighted UniFrac") +
scale_color_brewer(palette="Dark2")
plot_ordination(ps.lt.ts.nosurgery, LT_nosurgery.uni_ord, type = "samples", color = "Intranasal_Treatment") + theme_bw(base_size = 14) +
scale_color_brewer(palette="Dark2") +
theme(legend.position = "top") -> p.lt.nosurgery.uni
p.lt.nosurgery.uni
plot_ordination(ps.lt.ts.nosurgery, LT_nosurgery.weighted_uni_ord, type = "samples", color = "Intranasal_Treatment") +
theme_bw(base_size = 14) +
scale_color_brewer(palette="Dark2") +
theme(legend.position = "top") -> p.lt.nosurgery.wuni
p.lt.nosurgery.wuni
plot_grid(p.lt.nosurgery.uni, p.lt.nosurgery.wuni,
ncol = 2, labels = "AUTO", label_size = 20) -> p.lt.nosurgery
p.lt.nosurgery
ragg::agg_jpeg(file.path(fig.dir, "fig_lt_nosurgery_beta.jpeg"), width = 10, height = 5, units = "in", res = 600, scaling = 1.1)
p.lt.nosurgery
dev.off()## png
## 25.2 Did Surgery result in lung microbiome changes?
Only look at at HFD at week 13.
ps.LT.asv.ts %>% subset_samples(Diet == "High_Fat_Diet") %>% subset_samples(Week == "13") -> ps.lt.ts.hfd
# Calculate ditances
ps.lt.ts.hfd %>% phyloseq::distance(method = "unifrac") -> LT_hfd_uni
ps.lt.ts.hfd %>% phyloseq::distance(method = "wunifrac") -> LT_hfd_weighted_uni
# Ordinate distances
LT_hfd.uni_ord <- ordinate(ps.lt.ts.hfd, method = "PCoA", distance = "unifrac")
LT_hfd.weighted_uni_ord <- ordinate(ps.lt.ts.hfd, method = "PCoA", distance = "wunifrac")
# Extract dataframe
ps.lt.ts.hfd@sam_data %>% data.frame() -> LT_hfd.df
LT_hfd.df %>% dplyr::count(Surgery) %>% arrange(n)# Set seed for reproducibility
# UniFrac
set.seed(123)
adonis2(formula = LT_hfd_uni ~ Surgery, data = LT_hfd.df, permutations = 9999, parallel = getOption("mc.cores"))# Weighted UniFrac
set.seed(123)
adonis2(formula = LT_hfd_weighted_uni ~ Surgery, data = LT_hfd.df, permutations = 9999, parallel = getOption("mc.cores")) NS.
plot_ordination(ps.lt.ts.hfd, LT_hfd.uni_ord, type = "samples", color = "Surgery") + theme_bw(base_size = 14) +
scale_color_brewer(palette="Dark2") +
theme(legend.position = "top") -> p.lt.hfd.uni
p.lt.hfd.uni
plot_ordination(ps.lt.ts.hfd, LT_hfd.weighted_uni_ord, type = "samples", color = "Surgery") +
theme_bw(base_size = 14) +
scale_color_brewer(palette="Dark2") +
theme(legend.position = "top") -> p.lt.hfd.wuni
p.lt.hfd.wuni
plot_grid(p.lt.hfd.uni, p.lt.hfd.wuni,
ncol = 2, labels = "AUTO", label_size = 20) -> p.lt.hfd
p.lt.hfd
ragg::agg_jpeg(file.path(fig.dir, "fig_lt_hfd_beta.jpeg"), width = 10, height = 5, units = "in", res = 600, scaling = 1.1)
p.lt.hfd
dev.off()## png
## 26 Alpha Diversity: Fecal Samples
6.1 Calculate
ps.F.asv## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3304 taxa and 177 samples ]
## sample_data() Sample Data: [ 177 samples by 40 sample variables ]
## tax_table() Taxonomy Table: [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq() DNAStringSet: [ 3304 reference sequences ]ps.F.asv %>% estimate_richness(measures = c("Observed", "Shannon")) -> F_richness## Warning in estimate_richness(., measures = c("Observed", "Shannon")): The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
##
## We recommended that you find the un-trimmed data and retry.F_richness$Faith <- pd(otu_table(ps.F.asv), phy_tree(ps.F.asv), include.root = FALSE)$PD
F_richness %>% rownames_to_column(var = "Label") -> F_richness
gsub("X", "", F_richness$Label) -> F_richness$Label # remove X in front of sample names
meta.df %>%
filter(meta.df$Label %in% F_richness$Label) %>% # merge with metadata
select(c(Label, Genotype:Surgery, Kit:Feces_weight, Week, Cohort, Mouse, Notes)) %>%
right_join(F_richness, by = "Label") -> F_richness
F_richness %>% group_by(Mouse, Week) %>% filter(n()>1) -> F_richness.dup
F_richness.dup# calculate averages for duplicates
F_richness.dup %>%
dplyr::select(Label, Mouse, Week, Observed:Faith) %>% group_by(Mouse) %>%
dplyr::summarise(Observed = mean(Observed),
Shannon = mean(Shannon),
Faith = mean(Faith)) -> F_richness.dup.calc
# add week information for duplicates
F_richness.dup %>% dplyr::select(Mouse, Week) %>% unique() %>%
ungroup() %>%
right_join(F_richness.dup.calc, by = "Mouse") -> F_richness.dup.calc
# add metadata for duplicates
meta.df %>% filter(Sample_type == "Fecal") %>%
filter(Mouse %in% F_richness.dup.calc$Mouse) %>%
select(c(Genotype:Surgery, Week, Cohort, Mouse, Notes)) %>%
right_join(F_richness.dup.calc, by = c("Mouse", "Week")) %>%
unique() -> F_richness.dup.calc
F_richness.dup.calc# Clean up: remove duplicates from original dataframe, then add calculated values.
F_richness %>% filter(!Label %in% F_richness.dup$Label) %>%
bind_rows(F_richness.dup.calc) -> F_richnessF_richness %>% filter(is.na(Diet) == FALSE) -> F_richness
F_richness %>% pivot_longer(
cols = contains(c("Observed", "Shannon", "Faith")),
names_to = c("Measurement")
) -> F_richness_stat # dataframe that works well for ggplot2() and anova_test()
F_richness_stat$Mouse <- factor(F_richness_stat$Mouse)
F_richness_stat$Measurement <- factor(F_richness_stat$Measurement)
F_richness_stat$Cohort <- factor(F_richness_stat$Cohort)6.2 No Surgery group: effect of diet?
I want to see if high fat diet resulted in a significant decline in
alpha diversity, as has been reported before. To do this, I will compare
wk 0 and wk 10 alpha diversity for all mice
based on diet (to increase n). I also want to see if genotype
and sex had an effect.
F_richness %>% filter(Surgery == "None") -> F_richness.nosurgery
F_richness.nosurgery %>%
dplyr::count(Week, Diet, Genotype, Sex) %>% arrange(n)F_richness.nosurgery %>%
dplyr::count(Week, Diet, Genotype) %>% arrange(n)F_richness.nosurgery %>% distinct(Mouse) %>% nrow()## [1] 58F_richness.nosurgery %>% nrow()## [1] 117F_richness_stat %>% filter(Surgery == "None") -> F_richness.nosurgery.long
F_richness.nosurgery.long %>%
group_by(Measurement, Week, Diet, Genotype) %>%
shapiro_test(value) %>%
add_significance() -> F_richness.nosurgery.shapirowilk
F_richness.nosurgery.shapirowilkF_richness.nosurgery.shapirowilk %>% filter(p<0.05)All indistinguishable from normal distribution.
F_richness.nosurgery.long %>%
group_by(Measurement, Week, Diet, Genotype) %>%
identify_outliers(value)The feces weight, gDNA concentrations, and notes (or the lack thereof) for even extreme outliers look good.
F_richness.nosurgery.long %>%
group_by(Measurement) %>%
anova_test(dv = value,
wid = Mouse,
within = Week,
between = c(Diet, Genotype)) %>%
get_anova_table() %>%
as.data.frame() %>%
adjust_pvalue(method = "holm") %>%
add_significance() -> F.nosurgery.res
F.nosurgery.resF.nosurgery.res %>% filter(p.adj < 0.05)Let’s graph what they look like:
# Change facet label names
alpha.labs <- c("Faith's PD", "Observed Richness", "Shannon")
names(alpha.labs) <- c("Faith", "Observed", "Shannon")
ggplot(F_richness.nosurgery.long, aes(x = as.factor(Week), y = value, color = Diet))+
geom_boxplot(position=position_dodge(width = 0.8)) + geom_point(position=position_dodge(width = 0.8)) +
facet_wrap(~Measurement, scales = "free", labeller = labeller(Measurement = alpha.labs)) +
labs(x="Week", y="Alpha Diversity Metric Value", fill="Diet") + expand_limits(y = 0) + theme_bw() +
theme(legend.position = "top") 
# plot mean +/- standard error
fun_se <- function(x){
c(std_error = sd(x)/sqrt(length(x)), mean = mean(x))
}
# calculate mean and se for each week, measurement, and diet
F_richness.nosurgery.long %>%
doBy::summary_by(value ~ Measurement + Week + Diet, FUN = fun_se) -> F_richness.nosurgery.graph
ggplot(F_richness.nosurgery.graph, aes(x = Week, y = value.mean, color = Diet))+
geom_line(aes(group=Diet)) +
geom_errorbar(aes(ymin=value.mean-value.std_error, ymax=value.mean+value.std_error, width=1)) +
theme_bw(base_size = 14) + scale_fill_brewer(palette="RdBu") +
facet_wrap(~Measurement, scales = "free_y",
labeller = labeller(Measurement = alpha.labs)) +
labs(x="Week", y="Alpha Diversity Metric Value", fill="Diet") + expand_limits(y = 0) +
theme(legend.position = "top") -> fig.alpha.nosurgery
fig.alpha.nosurgery
ragg::agg_jpeg(file.path(fig.dir, "fig.alpha.nosurgery.jpeg"), width = 6, height = 5, units = "in", res = 600)
fig.alpha.nosurgery
dev.off()## png
## 26.3 Effect of surgery?
F_richness %>%
filter(Diet == "High_Fat_Diet" & Week == 13) -> f.richness.hfd
f.richness.hfd %>% dplyr::count(Surgery) %>% arrange(n)F_richness %>%
filter(Diet == "High_Fat_Diet" & Week == 13) %>% distinct(Mouse) %>% nrow()## [1] 45f.richness.hfd %>% dplyr::count(Surgery, Genotype) %>% arrange(n)f.richness.hfd %>% dplyr::count(Surgery, Genotype, Sex) %>% arrange(n)F_richness_stat %>%
filter(Diet == "High_Fat_Diet" & Week == 13) -> f.richness.hfdf.richness.hfd %>%
group_by(Measurement, Surgery, Genotype) %>%
shapiro_test(value) %>%
add_significance()f.richness.hfd %>%
group_by(Measurement, Surgery, Genotype) %>%
identify_outliers(value)f.richness.hfd %>% filter(Genotype == "KO") %>%
ggqqplot(x = "value") + facet_grid(vars(Measurement), vars(Surgery), scales = "free")
The QQ plots look quite good, except for the few outliers.
f.richness.hfd %>%
group_by(Measurement) %>%
anova_test(dv = value,
wid = Mouse,
between = c(Surgery, Genotype, Sex)) %>%
get_anova_table() %>%
as.data.frame() %>%
adjust_pvalue(method = "holm") %>%
add_significance() -> F.hfd.res
F.hfd.resAll ns. Let’s graph:
f.richness.hfd %>%
ggplot(aes(x=Surgery, y=value, color=Surgery)) +
geom_boxplot() +
geom_point() +
scale_fill_brewer(palette="RdBu") +
labs( y="Alpha Diversity Metric value", color="Surgery", x = "Surgery") +
facet_wrap(vars(Measurement), scales = "free", labeller = labeller(Measurement = alpha.labs)) +
theme_bw(base_size = 14) + expand_limits(y = 0) + theme(legend.position = "top") -> fig.alpha.hfd
fig.alpha.hfd
ragg::agg_jpeg(file.path(fig.dir, "fig_f_alpha_hfd.jpeg"), width = 7, height = 5, units = "in", res = 600)
fig.alpha.hfd
dev.off()## png
## 27 Beta Diversity with fecal samples: Diet
ps.F.asv.ts %>% subset_samples(Surgery == "None") -> F_Diet
F_Diet %>%
phyloseq::distance(method = "unifrac") -> F_uni
F_Diet %>%
phyloseq::distance(method = "wunifrac") -> F_weighted_uni
F_Diet@sam_data %>% data.frame() -> F.df
F_uni_ord <- ordinate(F_Diet, method = "PCoA", distance = "unifrac")
F_weighted_uni_ord <- ordinate(F_Diet, method = "PCoA", distance = "wunifrac")F.df %>% dplyr::count(Diet, Week) %>% arrange(n)F.df %>% dplyr::count(Diet, Genotype, Week) %>% arrange(n)F.df %>% nrow()## [1] 120F.df %>% distinct(Mouse) %>% nrow()## [1] 58# Set seed for reproducibility
set.seed(123)
adonis2(formula = F_uni ~ Diet * Genotype * Week, data = F.df, permutations = 9999, parallel = getOption("mc.cores")) -> F_diet_uni.res
F_diet_uni.resset.seed(123)
adonis2(formula = F_weighted_uni ~ Diet * Genotype * Week, data = F.df, permutations = 9999, parallel = getOption("mc.cores")) -> F_diet_weighted_uni.res
F_diet_weighted_uni.resAs expected, diet is highly significant. Genotype is significant but not any of its interaction terms.
F.df %>% dplyr::count(Genotype, Diet) %>% arrange(n)plot_ordination(F_Diet, F_uni_ord, type = "samples", color = "Diet") + theme_bw(base_size = 14) +
scale_fill_brewer(palette="RdBu") + theme(legend.position = "top") + labs(title = "UniFrac distance") +
theme(plot.title = element_text(hjust = 0.5)) +
facet_wrap(~Week) + stat_ellipse() -> fig.diet.uni
plot_ordination(F_Diet, F_weighted_uni_ord, type = "samples", color = "Diet") +
theme_bw(base_size = 14) +
scale_fill_brewer(palette="RdBu") + theme(legend.position = "top", panel.spacing.x = unit(10, "mm")) +
labs(title = "Weighted UniFrac distance") + theme(plot.title = element_text(hjust = 0.5)) +
facet_wrap(~Week) + stat_ellipse() -> fig.diet.w.uni
fig.diet.uni 
fig.diet.w.uni
plot_ordination(F_Diet, F_uni_ord, type = "samples", color = "Genotype") + theme_bw(base_size = 14) + labs(title = "UniFrac distance") +
scale_fill_brewer(palette="RdBu") + theme(legend.position = "top", plot.title = element_text(hjust = 0.5)) + stat_ellipse() -> fig.geno.uni
plot_ordination(F_Diet, F_weighted_uni_ord, type = "samples", color = "Genotype") + labs(title = "Weighted UniFrac distance") +
theme_bw(base_size = 14) +
scale_fill_brewer(palette="Dark2") + theme(legend.position = "top", panel.spacing.x = unit(8, "mm"), plot.title = element_text(hjust = 0.5)) + stat_ellipse() -> fig.geno.w.uni
# UniFrac: diet*genotype was significant
plot_ordination(F_Diet, F_uni_ord, type = "samples", color = "Genotype") + theme_bw(base_size = 14) +
scale_color_brewer(palette="Dark2") + facet_wrap(~Diet) + labs(title = "UniFrac distance") +
theme(legend.position = "top", plot.title = element_text(hjust = 0.5)) + stat_ellipse() -> fig.diet.geno.uni
fig.geno.uni
fig.diet.geno.uni
fig.geno.w.uni
plot_ordination(F_Diet, F_uni_ord, type = "samples", color = "Genotype") + theme_bw(base_size = 14) +
scale_fill_brewer(palette="RdBu") +
theme(legend.position = "top", panel.spacing.x = unit(5, "mm")) +
facet_grid(vars(Diet),vars(Week)) + stat_ellipse() -> fig.diet.geno.week.uni
fig.diet.geno.week.uni
plot_grid(
fig.alpha.nosurgery,fig.diet.uni,
fig.diet.w.uni, fig.diet.geno.uni, labels = "AUTO", ncol = 2, label_size = 20
) -> fig.diet.summary
fig.diet.summary
ragg::agg_jpeg(file.path(fig.dir, "fig_diet_summary.jpeg"), width = 18, height = 12, units = "in", res = 600)
fig.diet.summary
dev.off()## png
## 2plot_grid(
fig.geno.uni, fig.geno.w.uni, labels = "AUTO", ncol = 2, label_size = 20
) -> fig.nosurgery.geno.summary
fig.nosurgery.geno.summary
ragg::agg_jpeg(file.path(fig.dir, "fig_geno_summary.jpeg"), width = 12, height = 6, units = "in", res = 600)
fig.nosurgery.geno.summary
dev.off()## png
## 27.1 Fecal beta diversity: HFD
HFD: Effect of surgery; focus on wk 13 only.
ps.F.asv.ts %>% subset_samples(Diet == "High_Fat_Diet" & Week == 13) -> F_HFD
F_HFD## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3304 taxa and 45 samples ]
## sample_data() Sample Data: [ 45 samples by 40 sample variables ]
## tax_table() Taxonomy Table: [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq() DNAStringSet: [ 3304 reference sequences ]F_HFD %>%
phyloseq::distance(method = "unifrac") -> F_uni.HFD
F_HFD %>%
phyloseq::distance(method = "wunifrac") -> F_weighted_uni.HFD
F_HFD@sam_data %>% data.frame() -> F_HFD.df
F_uni_ord.HFD <- ordinate(F_HFD, method = "PCoA", distance = "unifrac")
F_weighted_uni_ord.HFD <- ordinate(F_HFD, method = "PCoA", distance = "wunifrac")F_HFD.df %>% dplyr::count(Surgery) %>% arrange(n)F_HFD.df %>% dplyr::count(Genotype) %>% arrange(n)F_HFD.df %>% dplyr::count(Surgery, Genotype) %>% arrange(n)F_HFD.df %>% nrow()## [1] 45# Set seed for reproducibility
set.seed(123)
adonis2(formula = F_uni.HFD ~ Surgery * Genotype, data = F_HFD.df, permutations = 9999, parallel = getOption("mc.cores")) set.seed(123)
adonis2(formula = F_weighted_uni.HFD ~ Surgery * Genotype, data = F_HFD.df, permutations = 9999, parallel = getOption("mc.cores")) Genotype is significant, surgery is not.
plot_ordination(F_Diet, F_weighted_uni_ord, type = “samples”, color = “Genotype”) + labs(title = “Weighted UniFrac distance”) + theme_bw(base_size = 14) + scale_fill_brewer(palette=“RdBu”) + theme(legend.position = “top”, panel.spacing.x = unit(8, “mm”), plot.title = element_text(hjust = 0.5)) + stat_ellipse() -> fig.geno.w.uni
# stat_ellipse for the ones that are *
plot_ordination(F_HFD, F_uni_ord.HFD, type = "samples", color = "Genotype") + theme_bw(base_size = 14) +
labs(title = "UniFrac distance") +
theme(legend.position = "top", plot.title = element_text(hjust = 0.5)) + stat_ellipse() -> fig.hfd.uni.geno
fig.hfd.uni.geno
plot_ordination(F_HFD, F_weighted_uni_ord.HFD, type = "samples", color = "Genotype") + theme_bw(base_size = 14) +
labs(title = "Weighted UniFrac distance") +
theme(legend.position = "top", plot.title = element_text(hjust = 0.5)) +
stat_ellipse() -> fig.hfd.w.uni.geno
fig.hfd.w.uni.geno
# Graphing out of curiosity
plot_ordination(F_HFD, F_uni_ord.HFD, type = "samples", color = "Surgery") + theme_bw(base_size = 14) +
scale_color_brewer(palette="Dark2") + theme(legend.position = "top") + stat_ellipse()
plot_ordination(F_HFD, F_weighted_uni_ord.HFD, type = "samples", color = "Surgery") + theme_bw(base_size = 14) +
scale_color_brewer(palette="Dark2") + theme(legend.position = "top") + stat_ellipse()
plot_grid(
fig.hfd.uni.geno,
fig.hfd.w.uni.geno,
labels = "AUTO", ncol = 2, label_size = 20
) -> fig.hfd.summary.beta
fig.hfd.summary.beta
ragg::agg_jpeg(file.path(fig.dir, "fig_hfd_beta.jpeg"), width = 12, height = 5, units = "in", res = 600)
fig.hfd.summary.beta
dev.off()## png
## 28 Reproducibility
sessionInfo()## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] ggtext_0.1.2 cowplot_1.1.1
## [3] ragg_1.2.6 doBy_4.6.19
## [5] DESeq2_1.40.2 SummarizedExperiment_1.30.2
## [7] Biobase_2.60.0 MatrixGenerics_1.12.3
## [9] matrixStats_1.1.0 GenomicRanges_1.52.1
## [11] GenomeInfoDb_1.36.4 IRanges_2.34.1
## [13] S4Vectors_0.38.2 BiocGenerics_0.46.0
## [15] phytools_1.9-16 maps_3.4.1
## [17] picante_1.8.2 nlme_3.1-163
## [19] ape_5.7-1 plotly_4.10.3
## [21] DT_0.29 ggrepel_0.9.4
## [23] microbiome_1.22.0 colorBlindness_0.1.9
## [25] RColorBrewer_1.1-3 pals_1.8
## [27] microViz_0.11.0 ggpubr_0.6.0
## [29] rstatix_0.7.2 vegan_2.6-4
## [31] lattice_0.22-5 permute_0.9-7
## [33] lubridate_1.9.3 forcats_1.0.0
## [35] stringr_1.5.0 dplyr_1.1.3
## [37] purrr_1.0.2 readr_2.1.4
## [39] tidyr_1.3.0 tibble_3.2.1
## [41] ggplot2_3.4.4 tidyverse_2.0.0
## [43] fs_1.6.3 phyloseq_1.44.0
## [45] magrittr_2.0.3
##
## loaded via a namespace (and not attached):
## [1] rstudioapi_0.15.0 jsonlite_1.8.7 farver_2.1.1
## [4] rmarkdown_2.25 zlibbioc_1.46.0 vctrs_0.6.4
## [7] multtest_2.56.0 RCurl_1.98-1.13 S4Arrays_1.2.1
## [10] htmltools_0.5.6.1 plotrix_3.8-2 broom_1.0.5
## [13] Rhdf5lib_1.22.1 rhdf5_2.44.0 gridGraphics_0.5-1
## [16] sass_0.4.7 bslib_0.5.1 htmlwidgets_1.6.2
## [19] plyr_1.8.9 cachem_1.0.8 igraph_1.5.1
## [22] lifecycle_1.0.3 iterators_1.0.14 pkgconfig_2.0.3
## [25] Matrix_1.6-1.1 R6_2.5.1 fastmap_1.1.1
## [28] GenomeInfoDbData_1.2.10 digest_0.6.33 numDeriv_2016.8-1.1
## [31] colorspace_2.1-0 crosstalk_1.2.0 textshaping_0.3.7
## [34] labeling_0.4.3 clusterGeneration_1.3.8 fansi_1.0.5
## [37] timechange_0.2.0 httr_1.4.7 abind_1.4-5
## [40] mgcv_1.9-0 compiler_4.3.1 microbenchmark_1.4.10
## [43] bit64_4.0.5 withr_2.5.1 doParallel_1.0.17
## [46] backports_1.4.1 BiocParallel_1.34.2 optimParallel_1.0-2
## [49] carData_3.0-5 ggsignif_0.6.4 MASS_7.3-60
## [52] DelayedArray_0.26.7 scatterplot3d_0.3-44 biomformat_1.28.0
## [55] tools_4.3.1 glue_1.6.2 quadprog_1.5-8
## [58] rhdf5filters_1.12.1 gridtext_0.1.5 grid_4.3.1
## [61] Rtsne_0.16 cluster_2.1.4 reshape2_1.4.4
## [64] ade4_1.7-22 generics_0.1.3 gtable_0.3.4
## [67] tzdb_0.4.0 data.table_1.14.8 hms_1.1.3
## [70] xml2_1.3.5 Deriv_4.1.3 car_3.1-2
## [73] utf8_1.2.4 XVector_0.40.0 foreach_1.5.2
## [76] pillar_1.9.0 vroom_1.6.4 splines_4.3.1
## [79] bit_4.0.5 survival_3.5-7 tidyselect_1.2.0
## [82] locfit_1.5-9.8 Biostrings_2.68.1 knitr_1.44
## [85] xfun_0.40 expm_0.999-7 stringi_1.7.12
## [88] lazyeval_0.2.2 yaml_2.3.7 evaluate_0.22
## [91] codetools_0.2-19 cli_3.6.1 systemfonts_1.0.5
## [94] munsell_0.5.0 jquerylib_0.1.4 dichromat_2.0-0.1
## [97] Rcpp_1.0.11 mapproj_1.2.11 coda_0.19-4
## [100] ellipsis_0.3.2 bitops_1.0-7 phangorn_2.11.1
## [103] viridisLite_0.4.2 scales_1.2.1 crayon_1.5.2
## [106] combinat_0.0-8 rlang_1.1.1 fastmatch_1.1-4
## [109] mnormt_2.1.1